Analytical realization of complex thermal meta-devices

Fourier’s law dictates that heat flows from warm to cold. Nevertheless, devices can be tailored to cloak obstacles or even reverse the heat flow. Mathematical transformation yields closed-form equations for graded, highly anisotropic thermal metamaterial distributions needed for obtaining such functionalities. For simple geometries, devices can be realized by regular conductor distributions; however, for complex geometries, physical realizations have so far been challenging, and sub-optimal solutions have been obtained by expensive numerical approaches. Here we suggest a straightforward and highly efficient analytical de-homogenization approach that uses optimal multi-rank laminates to provide closed-form solutions for any imaginable thermal manipulation device. We create thermal cloaks, rotators, and concentrators in complex domains with close-to-optimal performance and esthetic elegance. The devices are fabricated using metal 3D printing, and their omnidirectional thermal functionalities are investigated numerically and validated experimentally. The analytical approach enables next-generation free-form thermal meta-devices with efficient synthesis, near-optimal performance, and concise patterns.

1.In Figure 1A, the color of "Thermal conductivity distribution" is confused.According to transformation thermotics, the thermal conductivity is anisotropic.So what does the color mean?
We thank the reviewer for pointing this out.We have modified Figure 1A to reflect the anisotropy by illustrating the distribution of the conductivity tensor components as follows.
2. The process from "Thermal conductivity distribution" to "Analytical inverse homogenization" (Figure 1A) has not been described clearly in the manuscript.The thermal conductivity tensor relationship between these two parts is not clearly explained.
We thank the reviewer for pointing this out.We have improved the description as follows.
"The eigenvalues are used for determining the volume fractions of Material A at the two scales of the rank-2 laminate ( 1 and  2 ) through closed-form inverse homogenization, which is stated as finding  1 and  2 such that the resulting laminate's homogenized principal conductivities are equal to the eigenvalues (see SI 3.1 and SI 3.2 for details)." 3. In Figure 1B, the variate of the equation in "Transformation thermotics" part should be r and θ.where physical values are  1 ,  2 ∈ [0,1] .In (11),  2, and  2, are set equal to the two eigenvalues of ." 5. It is still not clear how the final structure is designed with specific feature size and the other information that is input to the wave function.The authors should explain in detail how to get the final structure with pseudo-conformal mapping and feature size.Please take thermal cloak as an example to give a specific derivation.
We thank the reviewer for the suggestion.To comprehensively illustrate the process of generating the structures, we have added a section SI 5.3 and a figure (Figure 5) to the Supporting Information as follows.
"After obtaining the local microstructural parameters  1 and  2 (in 3.3) and the mapping  1 and  2 (in 5.1 and 5.2) as four scalar fields on the meta-device domain, we can generate the corresponding global structure through a wave function representation [5].Specifically, we need two wave functions to represent the two groups of orthogonal laminates in the global structure, respectively.The wave functions are where the feature size parameter   is user-defined and needs to satisfy the mild requirement discussed in 5.2.
To illustrate the procedure, we use the cloak as an example as shown in Figure 5.Given the four (discrete) fields of  1 ,  2 ,  1 , and  2 , we obtain the wave functions  1 and  2 using (20) with values plotted in the third column of Figure 5. Based on [5], the Material 1 (steel) parts for the first and second laminates are defined by the regions  1 ≤ 0 and  2 ≤ 0, respectively, which are shown in the fourth column of Figure 5.Note that relative widths and orientations of the laminates are consistent with the distributions of   and   by construction.Finally, the metamaterial structure is obtained as the union of the two laminates as shown in the rightmost plot of Figure 5, and the union set can be expressed as   ≤ 0 with   ≔ min{ 1 ,  2 } .We want to emphasize that the process of generating the structure, i.e., computing the wave functions and taking the union, requires negligible computational cost and is free of post-processing for structural connectivity, which is a major advantage over state-ofthe-art computation-based approaches." 6.In figure 5A, why there is no isotherm in thermal cloak but the other two include it?
We thank the reviewer for pointing this out.We have added the isotherms to Figure 5A as follows.
We have also modified the description about the cloak experiment as follows.
"For the cloak, the temperature gradient in the background is uniform and largely unperturbed while the temperature gradient inside the core is very small as demonstrated by the isotherms."We thank the reviewers for suggesting these relevant studies.They have been cited in the first two paragraphs of the revised Introduction as follows.
"Metamaterials exhibit extreme properties not found in nature and enable exotic functionalities once conceivable only in fiction, such as optically cloaking an object from its environment [1,2,3].Originating in theoretical electromagnetism [1], the idea of transformation-based omnidirectional field manipulation through strategically designed material distribution quickly gained enormous attention and spread to multiple disciplines.Most prominently, cloaking of physical fields [4,5] are theoretically investigated and experimentally reproduced in the discipline of electromagnetism [4,2,3,6,7,8,9] elasticity [10,11,12], acoustics [13,14,15,16,17], and thermotics [18,19,20,21,22].In the latter, exotic functions such as heat cloaking [23,24,25,26,27], rotating [28], concentrating [29], camouflaging [30,31], illusion [32], and encrypting [33] have been designed and reproduced experimentally.In addition to transformation-based theories, direct use of topology optimization [34,35,36] and data-driven methods for designing global structures [37,38,39,40] or local microstructures [41] can produce the desired functionalities for specific targeted load cases and directions and are thus non-omnidirectional.We note that this class of non-omnidirectional metamaterials does not employ transformation thermotics, and the meta-device's performance significantly deteriorates when the direction of the applied heat is changed.The scope of this study is the omnidirectional class of thermal metamaterial, where furthermore the conductivity distribution is obtained analytically and free of iterative procedures at a fraction of the cost of abovementioned numerical approaches." Reviewer #2 (Remarks to the Author): In this paper, microstructures are designed by evaluation using finite element analysis and optimization using Newton's method to achieve the thermal properties required by transformation thermotics for realizing thermal cloak and related metadevices.Although it is interesting overall for this reviewer, it is difficult to recommend that the presented study is innovative enough to be published in Nature Communications.In addition, it is very difficult to understand the position of this research in optimization research on thermal cloak (and its applications), and the paper simply claims to be the state of the art.The reasons why this paper is NOT sufficient for publication in Nature Communications are described below.
1. IMPACT on METHOD: Despite many papers have been published on thermal cloaking (and its applications) designed by numerical approaches including topology optimization, but there are almost no descriptions of them.The authors should comprehensively cite and organize previous studies using the density method, level set method, and homogenization method (including papers using dehomogenization), and explain the location and advantages of this research.However, no descriptions and paragraph on the previous works are written in the introduction.
We thank the reviewer for the comments.We would like to clarify that in the original draft, we have cited and discussed in detail three papers [42,43,26] that use density-based topology optimization for designing thermal metamaterials in the Introduction.To address the concern, we have added comments and citations on level-set topology optimization approaches to thermal metamaterials.We have also included an additional paper using density-based topology optimization and discussed its features in the third paragraph of the Introduction.
As discussed in the Introduction of our manuscript, our study fundamentally departs from the state-of-the-art of both the non-omnidirectional (elaborated in this response comment) and omnidirectional thermal metamaterial (elaborated in the response comment #2 below).We would like to clarify that most topology optimization-based meta-devices belong to a fundamentally different non-omnidirectional class of meta-devices with load case-and direction-sensitive performance.That class of non-omnidirectional meta-devices does not use transformation thermotics, and performance significantly deteriorates when the direction of the applied heat is changed.That class of metamaterial is not the focus of this work and hence is not the main focus of our Introduction discussion (although also cited).We have clarified this point in the added paragraph.
"In addition to transformation-based theories, direct use of topology optimization [34,35,36] and data-driven methods for designing global structures [37,38,39,40] or local microstructures [41] can produce the desired functionalities for specific targeted load cases and directions and are thus non-omnidirectional.We note that this class of non-omnidirectional metamaterials does not employ transformation thermotics, and the meta-device's performance significantly deteriorates when the direction of the applied heat is changed.The scope of this study is the omnidirectional class of thermal metamaterial, where furthermore the conductivity distribution is obtained analytically and free of iterative procedures at a fraction of the cost of above-mentioned numerical approaches." "If not solvable by simple intuition, the state-of-the-art strategies to address the dehomogenization challenge is to divide the irregular domains into many square unit cells and use computational morphogenesis approaches (such as topology optimization) with numerical homogenization to inversely design the locally graded unit cell microstructures [42,43,26,44]." 2. Also, this study is very similar to the previous study [24], and the introduction need to include how this study is innovative compared to the study ref.[24] that used topology optimization using the density method for microstructure design.Of course, there are some differences from the references.However, when asked whether this has enough impact to be published in Nature Communications, this reviewer thinks it is NOT.
This study is fundamentally different from, and much advantageous over [24] (new reference ordering is [42]) in the following aspects: the local microstructure generation, global structure realization, structure connectivity, and thermal function performance, as evidenced by the following: o Straightforward analytical realization v.s.computational inversely designed microstructures: The rank-2 microstructures in our structures are obtained analytically with closed-form solutions.There is no design aspect involved in the entire process.The structures in [24] are obtained by topology optimization inverse design involving numerous heavy numerical solutions of PDEs and iterative optimization updates.
o Optimal v.s.sub-optimal performance: In our study, the closed-form solution of rank-2 laminates ensures that any given (physically achievable) target conductivity is achieved exactly, thereby ensuring the close-to-optimal performance.By contrast, the microstructures generated by topology optimization in [24] are sub-optimal and may not exactly achieve a given target due to the severe non-convexity of the inverse design problem as well as connectivity issues between cells.
o Guaranteed and post-processing-free structural connectivity v.s.disconnected and postprocessing-requiring structures: Our approach of using a pseudo-conformal mapping to generate the global structure naturally ensures connectivity among the members and is free of post-processing, whereas the global structure connectivity in [24] cannot be satisfied without heavy post-processing that in turn perturbs the local homogenized properties.
o Negligible v.s.extremely high computational cost: The main computational cost of our approach for a high-resolution structure is less than 30 seconds on a desktop workstation, whereas the approach in [24] requires orders-of-magnitude higher computation cost due to the need for solving many inverse design problems, each of which in turn requires several-hundred numerical solves of the PDE.
o Elegant and aesthetically concise v.s.hard-to-interpret and overly complex patterns: Our global structural patterns (with close-to-optimal performance) are elegant as the member orientations are aligned with the principal directions of the theoretical conductivity and member size proportional to the magnitude of eigenvalues of the conductivity.Further, our global structure is aesthetically concise because each local microstructure (rank-2 laminates) contains only two members.By contrast, the global structures in [24] have many more complex structural branches with highly irregular sizes, shapes, and orientations, which prohibit direct perception of the magnitude and principal directions of the conductivity and complicates the fabrication.
o On a higher and more general level, our study demonstrates an inspirational idea: complex and anisotropic material property distributions can be analytically and optimally realized by extremely simple graded structures.By contrast, study [24] appears to have demonstrated the reverse: realization of complex material properties requires highly sophisticated techniques and microstructures.The new inspiration from our study is translatable and impactful to many other engineering disciplines.
The above facts demonstrate that the study is fundamentally different from [24] and is not "very similar to the previous study" as suggested in the comment.All the above arguments were already included in our original manuscript, however, in light of the reviewer's comments, we have carefully gone through the manuscript and edited and sharpened, where appropriate.
3. PERFORMANCE: Despite the very high degree of design freedom by designing the microstructure, the RTD of the designed structure is higher (worsen) than the equivalent value in some previous studies and it has not been improved.We would like to clarify that the scope of this paper is the transformation-based, omnidirectional class of thermal metamaterial with performance insensitive to applied heat direction and material of the core (core can be conductive or insulating).Thus, the above comparison to the state-of-the-art in the omnidirectional class is valid.Our scope is not the direction-sensitive class used in many numerical works.For that class of metamaterial, the performance can be very good, if not perfect (= 0), exactly in the targeted direction and with an insulating core.However, change of heat direction or core material significantly deteriorates the performance to "inferior" as described in [37], with RTD measures increased to 15% -80% [37].Therefore, we do not (and should not) compare our performance with that fundamentally different and simpler class of meta-devices.
To avoid confusion about the metadevice class, we have added the description below to the first paragraph of the Introduction.
"In addition to transformation-based theories, direct use of topology optimization [34,35,36] and data-driven methods for designing global structures [37,38,39,40] or local microstructures [41] can produce the desired functionalities for specific targeted load cases and directions and are thus non-omnidirectional.We note that this class of non-omnidirectional metamaterials does not employ transformation thermotics, and the meta-device's performance significantly deteriorates when the direction of the applied heat is changed.The scope of this study is the omnidirectional class of thermal metamaterial, where furthermore the conductivity distribution is obtained analytically and free of iterative procedures at a fraction of the cost of above-mentioned numerical approaches." 4. QUALITY: (??) remains just before the SI equation (15) and this shows the quality of the manuscript.
We apologize for this editing error.We have corrected it and thoroughly rechecked the revised manuscript.
Reviewer #3 (Remarks to the Author): The so called "transformation physics" leads to designs, using composites at the microscale, that lead to sometimes surprising conclusions about what is possible.However, the designs are typically far from optimal and sometimes utilize constituent materials with extreme properties.By contrast the authors show how one can use relatively simple designs, incorporating rank 2 laminates at the microscale, that have a continuous morphing of structure on the mesoscale, and they verify their near optimal performance experimentally.This is a major advance and, although currently limited to two-dimensions, there is no reason that the work cannot be extended to three dimensions.However, a major flaw in the paper is that it perpetuates misconceptions about original discoveries.
omnidirectional class of thermal metamaterial, where furthermore the conductivity distribution is obtained analytically and free of iterative procedures at a fraction of the cost of abovementioned numerical approaches." 2. On the second page it is stated that "the complete spectrum of composite's conductivity can indeed be fully and analytically achieved by simple rank-2 laminates".Nowhere here is started the all important fact that this is limited to two-dimensional composites (and that rank-3 laminates are needed in three dimensions).The references here should include those of Tartar and Murat (which can be found in books on the theory of composites).Otherwise it is a bit unfair, as Professor Kohn, on a trip to the Soviet Union, had communicated the results of Tartar and Murat to Lurie and Cherkaev who incorporated them in their work, without attribution (they claim they had proved the bounds in an alternative way but utilized the approach of Murat and Tartar as it was simpler).Lurie and Cherkaev had independently obtained the two dimensional bounds.
We thank the reviewer for alluding us to this historical insight.We have revised the statement to be rigorous and cited the suggested works [48,49] as well an additional [50] as shown below.
"Contrary to the indication in [42,43] that the anisotropic and heterogeneous conductivities are difficult to realize by layered structures, the complete spectrum of a two-dimensional (2D) composite's conductivity can indeed be fully and analytically achieved by simple rank-2 laminates [46,47,48,49,50]." 3. The authors talk about cloaking and concentrating fields between a hot plate and a cold plate.
For the case where the hot and cold plates are fixed, this was solved using rank one laminates by Gibiansky, Lurie, and Cherkaev and Gibiansky in the paper Zhurnal tekhnicheskoi fiziki 58 (1988) 67-74, english translation in Sov.Phys.Tech.Phys.33 :38-42 (1988).This can be made more clear if one considers the periodic extension of their solution.So that work needs to be cited.
We thank the reviewer's suggestion.We have added the citation [56] (the original version) to the result section Sec.2.2 of the revised manuscript as follows."Higher concentration rates require more extreme and anisotropic  values as manifested by the larger radial members and smaller circulating members.As a side note, related problems for uni-directional concentration have been solved and require only rank-1 laminates for their realization [56]." 4. A major deficiency is that I do not see the equations of thermal conduction used by the authors stated anywhere.Presumably the authors are not incorporating Joule heating which is a non-linear term, but without some statement of the equations it is hard to know what they are including.
We have added the equation with the description of the setup in the first paragraph of Results as follows.
"This study focuses on steady-state heat conduction in a 2D medium without heat sources or sinks, which is mathematically described by ∇ ⋅ (∇) = 0 with  being the temperature field and  the anisotropic 2x2 temperature-independent thermal conductivity tensor of the medium."Additionally, we have added a section (SI 6) in the Supplementary Information about the thermal conduction boundary value problem and the finite element solution setup as follows.
"After obtaining the design of the metamaterial structure, we post-evaluate its heat conduction response using FEM.The steady-state 2D heat conduction problem is where Ω is the domain, ∂Ω , and ∂Ω , are the Dirichlet boundaries of the hot and cold ends, respectively, and ∂Ω  are the adiabatic Neumann boundaries, such that �∂Ω , ∪ ∂Ω , � ∪ ∂Ω  = ∂Ω and �∂Ω , ∪ ∂Ω , � ∩ ∂Ω  = ∅ ,  �  and  �  are the prescribed temperature for the hot and cold ends, respectively.Note that we focus on the case with no internal heat sources or sinks.
The boundary value problem is solved using an in-house FEM program where the four-node quadrilateral element with four Gauss points is adopted.The  is assumed to be piece-wise constant and associated with the elements.Specifically,  =  Steel  for the steel parts of the structure,  =  PDMS  for the PDMS parts, and  =  0  for the background.The squareshaped domain is discretized by a regular 2000x2000 mesh to allow for sufficient resolution of the microstructure."

General comment:
In the revision process, we have carefully reread the entire manuscript and made changes throughout that we believe have sharpened our messages.
The authors are grateful for the reviewers' recognition and additional insightful comments.We have addressed the comments and revised accordingly.The textual changes are highlighted in red color in the revised manuscript and the added citations are highlighted yellow.Below we address each comment in detail.

Reviewer #2 (Remarks to the Author):
I read the Response and reconsider the importance of the presented originality in the manuscript.
Major Comments: 1. IMPACT & PERFORMANCE: Indeed, as stated in authors' reply to my comment 2, the method presented in this manuscript has several advantages compared to the Ref. [42].On the other hand, if my understand is correct, all the above advantages is achieved by sacrificing the degree of freedom in designing the global structure described by wave function representation in Eq. ( 2).Although it is true that Ref. [42] has a high computational cost and may not have sufficient performance compared to the results of this manuscript, it also has an extremely high degree of freedom in design, which is one of the original academic appeals of metamaterials.
The main premise is that a omni directional thermal cloak that can make object undetectable from heat flux for all directions can be realized with a very simple structure like a bi-layer cloak in Ref. [23], so it is not difficult to realize its function.Of course, by normalizing the equation, its performance is also independent of the magnitude of the heat flux.
In that sense, the proposed method is inferior to Ref. [42] in terms of design freedom, and inferior to bilayer cloaks in terms of structural simplicity (ease of fabrication) and performance.For numerical performance evaluation of the bilayer cloak, see the end of Ref. [37].
The method presented in this manuscript achieves the above advantages by sacrificing the design freedom, but are they worthy of publication in Nature Communications?
We thank the reviewer for raising these important points.We would like to clarify that this work did not sacrifice any design freedom, because the rank-2 laminate has the full design freedom in that it achieves all physically possible thermal conductivities for a 2D composite, as indicated in SI Figure 2 (attached below).This has been theoretically proved [48 -52].The geometric simplicity of rank-2 laminate is an advantage (rather than a disadvantage) over the complex microstructures, and it is not "inferior to Ref. [42] in terms of design freedom" nor has it "sacrificed design freedom".
To clarify this point, we have modified the comment on rank-2 laminates in the second paragraph of the Results section as follows: "With their simple microstructures, rank-2 laminates have been shown to cover the entire range of physically achievable thermal conductivities for 2D composites [52], and therefore, our approach does not sacrifice any design freedom." For the bi-layer cloak in [23], we are aware that it can achieve cloaking but only for very simple geometry, such as the circular/spherical domain and not for general domains.The scope and focus of this study have always been cloaks, rotators, and concentrators with general and complex geometries, which has been a major challenge in the field because of the required highly anisotropic and graded conductivities as mentioned in recent studies such as [26,44,45].Cloaks of general and complex geometries are thus not "not difficult to realize", and they are effectively and innovatively achieved in this study.Our study tackles a fundamentally different and more general scenario that cannot be realized by the bi-layer scheme [23].

Minor Comments
2. What is \rho used to express the geometry in Eq. ( 2) ?According to Ref. [53], \rho seems to be a level set function, but is this not the case?There might be no definition.
In the state-of-the-art paper [Y.Wang, W. Sha, M. Xiao, C. Qiu, L. Gao.Adv.Mater.2023, 2302387] on free-form omnidirectional thermal metamaterials, the Relative Temperature Difference (RTD) values of irregular-shape cloaks are higher than 13 % (the regular circular shape has 2.02%).By contrast, using the same RTD definition and with a similar level of shape complexity and irregularity, our corresponding omnidirectional cloaks achieve RTD values of 1.37% and 1.27% (as demonstrated in the figure below), which are much smaller than the current state-of-the-art irregular cloaks in [Y.Wang, W. Sha, M. Xiao, C. Qiu, L. Gao.Adv.Mater.2023, 2302387].We note that 0% RTD is the theoretically ideal situation, which cannot be attained by physical materials.By comparing the RTD values, we demonstrate a significant improvement over the state-of-the-art work on omnidirectional thermal metamaterials, and are not "worse than other papers" as claimed in Reviewer 2's comment about the RTD values.